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Abstract. We present a detailed analysis of the evolution of a simulated isolated disc galaxy. The simulation 
includes stars, gas, star formation and simple chemical yields. Stellar particles are split in two populations: the 
old one is present at the beginning of the simulation and is calibrated according to various ages and metallicities; 
the new population borns in the course of the simulation and inherits the metallicity of the gas particles. The results 
^ . have been calibrated in four wavebands with the spectro-photometric evolutionary model GISSEL2000 (Bruzual 

' & Chariot I1993II . Dust extinction has also been taken into account. A rest-frame morphological and bidimensional 

, photometric analysis has been performed on simulated images, with the same tools as for observations. The effects 

of the stellar bar formation and the linked star formation episode on the global properties of the galaxy (mass 
' and luminosity distribution, colours, isophotal radii) have been analysed. In particular, we have disentangled the 

, effects of stellar evolution from dynamic evolution to explain the cause of the isophotal radii variations. We show 

\l ■ that the dynamic properties (e.g. mass) of the area enclosed by any isophotal radius depends on the waveband 

and on the level of star formation activity. It is also shown that the bar isophotes remain thinner than mass 
isodensities a long time (> 0.7 Gyr) after the maximum of star formation rate. We show that bar ellipticity 
^ l' is very wavelength dependent as suggested by real observations. Effects of dust extinction on photometric and 

morphological measurements are systematically quantified. For instance, it is shown that, when the star formation 
rate is maximum, no more than 20% of the B band luminosity can escape from the bar region whereas, without 
^ dust extinction, bar B band luminosity accounts for 80% of the total B band luminosity. Moreover, the extinction 

, is not uniformly distributed inside the bar. 
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1. Introduction timescale (less than 1 Gyr). For instance, the gas can be 

responsible for the bar within bar phenomenon (Friedli 

Early N-body models of isolated disc galaxies, using only & Martinet ITTMIl or the dissolution of the bar and the 

collisionless particles, have clarified how several morpho- formation of a bulge by secular evolution (Pfenniger & 

logical features such as bars, spiral arm s (e g Hohl^978J Norman ITO Norman et al. lMll . Studies of the fuehng 

and boxy bulges (Combes & Sanders [Mil), appear in mechanisms at work in central regions (e.g. FriedU & Benz 

stellar discs and evolve. Soon afterwards, it was realised r^MI Shlosman & Nogiichi [Tm^ have been also studied 

that the role of the gas, the dissipative component, can- -v^jth such hybrid codes 
not be neglected. Collisionless N-body codes were thus 

coupled with hydrodynamic codes , either for cosmologi- However, stars and gas are not only bound together by 

cal purposes (Katz & Gunn 11 ^ , or for detailed stud- gravitation, but also by the process of star formation and 

ies of isolated (Friedli & Benz [1993|) or merging galaxies the energy and chemical feedback of supernova; and stellar 

(Mihos & Hernquist|l996j. For isolated disc galaxies, this ^inds. Since star formation is not well understood even in 

kind of codes has brought clear evidence that the pres- our own Galaxy, several recipes have been implemented 

ence of gas can deeply change the morphology on a short to mimic star formation in hydro-f stellar dynamic codes. 

These recipes have allowed to study for the first time the 

Send offprint requests to: Leo Michel-D ansae, e-mail: influence of bars on galactic abundance gradient (Friedli 

leo.michel-dansac@oamp.fr et al. lift!) 4)1 . the impact of the Stellar ultraviolet radiation 
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on star formation (Gerritsen & Icke I1997|l . the nuclear 
activity (Heller & Shlosman 1 1994(1 . the starbursts induced 
by mergers (e.g. Mihos & Hernquist '1996^1, the renewal 
of bars by gas accretion (Bournaud & Combes r2UU2|l . etc. 
They also allow to perform self-consistent simulations of 
the formation and evolution of galaxies (e.g. Katz 119921 
Steinmetz & Muller [TM5|l . 

Whatever the kind of code used for the study of secular 
evolution, the morphological analysis of the simulations fo- 
cuses on the properties of the mass distribution. To com- 
pare in details these models with observations it is implic- 
itly assumed that all particles have the same mass-to-light 
ratio taken to unity for convenience. However, the real 
galaxies are built of composite stellar populations, then 
of various mass-to-light ratios and metallicities. Thus, in 
order to make more straightforward comparisons of simu- 
lations with multi-wavelength observations, we have used 
stellar population synthesis models to photometrically cal- 
ibrate our self-consistent simulations which include stars, 
gas and star formation. We also take into account the ef- 
fect of dust extinction. 

This approach has been recently used to study the 
formation and the subsequent evolution of elliptical (e.g. 
Kawataf20Uli or disc galaxies in a A cold dark matter sce- 
nario (e.g. Westera et al.EnOl Abadi et al HUU^ . Our own 
approach slightly differs from previous ones since we do 
not simulate the formation of a disc galaxy in a primeval 
Universe. Indeed, we restrict ourselves to make a detailed 
study of the short-term evolution of an isolated disc galaxy 
already formed. All the particles are thus used to simulate 
the evolution of the galaxy. 

We describe in Sect. 2 the numerical model and our 
technique of photometric calibration. The global evolution 
(i.e. integrated properties) of the model is presented in 
Sect. 3, and the morphological evolution in Sect. 4. We 
summarise our findings in Sect. 5. 

2. The N-body simulation and photometric 
calibration 

2.1. Description of tfie code 

We used PMSPHSF, the N-body code developed in Geneva. 
It includes stars, gas and recipes to simulate star forma- 
tion. The broad outlines of the code are the following: 
the gravitational forces are computed with a particle- 
mesh method using a 3D polar grid with {Njf , N^, Nz) = 
(31, 32, 64) active cells (Pfenniger & Friedh 1993), the hy- 
drodynamics equations are solved using the SPH tech- 
nique (Friedli & Benz 11993(1 . Star formation recipes are 
those implemented by Friedli & Benz 1(1995(1 . Here we 
briefly summarise the main features. For a detailed de- 
scription of this modelling and a discussion on the influ- 
ence of the free parameters, we refer the reader to Friedli 

& Benz (nnns)- 

Star formation and feedback modelling is basically 
based on the instantaneous star formation approximation. 
It consists in the following steps : the identification of 



the regions of star formation, the conversion of a frac- 
tion of gas into stars and the computation of the amount 
of energy and heavy elements injected in the interstellar 
medium (i.e. the energy and chemical feedback). Because 
the latter step leads to the heating of gas, a simple treat- 
ment of radiative cooling is implemented into the energy 
equation. 

The first task is to identify the gaseous particles that 
will form stars. Friedli & Benz (1993) examined some pos- 
sible criteria. Not surprisingly, they found that the stan- 
dard Jeans instability criterion can be applied to spher- 
ical non-rotating gaseous systems, but that for rotating 
flat discs Toomre's instability criterion (Toomre 11964(1 
is a better indicator. Observational evidence also appear 
to support the use of this criterion as a good indicator 
for locating star formation at intermediate scale-lengths 
(Kennicutt IT^ |T990) . 

In any case, a particle i will be assumed to undergo a 
star formation episode if the following condition is verified: 



where is Toomre's parameter, Si is the local sound 
speed, Ki is the generalised epicyclic frequency, is the 
gas surface density. This condition is clearly Toomre's sta- 
bility criterion. The constant A equals unity in the case of 
an axisymmetric gaseous disc subject to radial instabili- 
ties (cf. Toomre Il964() . However, a value of A « 1.4 (with 
s = 6 km s~^) is derived from observations (Kennicutt 

Once a particle of mass subject to star formation 
has been identified, a fraction of its mass is converted into 
a stellar particle while total energy and momentum are 
conserved. The fraction of the gas particle actually trans- 
formed (i.e. the star formation efficiency factor e) is a free 
parameter. Note that within the framework of instanta- 
neous star formation, e is the product between the frac- 
tion of gas actually transformed into star 77 (in accord with 
the IMF) and the mass fraction 1 — R finally blocked in 
low mass stars or stellar remnants. Thus, we can write 
e = (1 — i?)?7. We re-use values of Friedh & Benz 1(1995(1 
which are R k. 0.45 (Maeder ITU^ and 0.02 < ry < 0.5. 
This translates into a possible range 0.01 < e < 0.25. For 
the present work, e = 0.1. 

In practice, a particle i subject to star formation is 
split into one star particle m*fj"'^"' and one gas particle 
^s,neio ^£ j.g(jy(,g(j niass. Mass conservation requires that 
these masses be given by 

*,new q q.new \ q 

TTij. =em", ml =(1 — e)m". 

At this stage, note that no Schmidt law (Schmidt 11959(1 
has been imposed; the star formation rate (SFR) is thus a 
consequence of the (hydro)dynamic instabilities (see Sect. 
5.2 of Friedli & Benz .1995 for a discussion on the rela- 
tion between the SFR and the gas surface density). The 
system can make several generations of stars from one 
gas particle and can thus quickly adjust itself to stabilise 
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any collapsing region. However, in order to avoid numeri- 
cal problems associated with the existence of gas particles 
differing widely in mass we restrict the allowable gas par- 
ticle mass range to a maximum of 20. For e = 0.1, this 
corresponds to a maximum of 29 generations of new stars 
from a single gas particle. Gas particles which reach the 
lower mass limit are transformed into a single coUisionless 
particle. 

Probably the most important effect of star formation 
on the global evolution of a galaxy is connected to the 
large amount of energy released in supernova explosions 
and stellar winds. In this modelling of feedback, the con- 
tribution of SNIa supernovae is neglected. Let Ef^f be 
the total energy produced by SNII explosions at particle 
z at a given timestep 



pSN 



where Eff^ - (1 - f)El^i, and i?,?^,, = are re- 

spectively the thermal and mechanical energy injected, 
R^'^ is the number of SN created par unit mass of stars 
formed, E^^ w lO^^ erg is the mean energy produced 
by one SNII, and / is a very poorly determined free pa- 
rameter between and 1. It can be seen as the transfer 
efficiency of mechanical energy by thermalization to the 
surrounding ISM. For the present work, we use / = 0.01 
and R^^ = 0.004. 

Depending on the model assumptions (Z, IMF, SFR), 
the asymptotic value of energy ratio R^ between hot 
winds and SNII contributions can vary from almost to 
about 2 (Leitherer et al. [1992). In particular, an increase 
of metals strongly increases R^ . So, the contribution of 
winds to the energy injection is simply assumed to be pro- 
portional to the contribution coming from SNe but with a 
linear metallicity dependence, i.e E^tot — R^ ^^ftot with 
i?^ = 2hZ. 

Energy, in both forms, is injected instantaneously and 
is spread over the ni neighbour particles j of particle i us- 
ing the SPH smoothing kernel as the normalised weighting 
function Wij . The fraction of mechanical energy given to 
particle j, Wij{Ef^^^^ + E^mec)^ is added by computing 
the corresponding velocity vector increment 



{Av 



l2W,,f{l + R^)EfZ r,-r, 



Similarly, the specific thermal energy E^^ increment re- 
ceived by particle j is given by 

^^_W,,{l~f)il + R^)Eff^, 



Au 



J 



1. 



Since the instantaneous star formation approxima- 
tion is used, instantaneous recycling is assumed as well. 
Therefore, during evolution, gas particles will see their 
abundances modified as star formation proceeds in either 
immediate neighbourhood, whereas star particles will be 
endowed at birth with a relative abundance pattern cor- 
responding to their parents gas cloud. Changes of star 



abundance due to their evolution are neglected. Formally, 
the abundance change of a gas particle over one timestep 
can be written as 



znt + At)^zut) 



where Wij is the weighting factor described previously, 
and yz is the net yield for metals (i.e. the mass of metals 
ejected by all stars per unit mass locked into stars). The 
net yield are taken from Maeder H1992L model C, Table 7). 
At birth, star particles are created from gas, thus they 
have the same chemical composition. 

The method of treating radiation losses is to use tabu- 
lated cooling functions from published data and to inter- 
polate between them. For the present work, we use those of 
Boehringer & Hcnslcr ( 1989 ) They cover the temperature 
range 10** K to 10* K. Therefore, we assume a lower cut-off 
of Tcut — 10^ K. Since we are not in a position to com- 
pute accurately the evolution of all elements, we restrict 
ourselves to cooling functions assuming solar abundance. 

2.2. Initial conditions 

An initial stellar population is setup to reproduce a disc 
galaxy with an already formed bulge. These particles form 
what we call hereafter the 'old population' as opposed to 
particles created during the evolution ('new population'). 

The initial stellar positions for 500 000 particles 
are drawn from a superposition of two axisymmetric 
Miyamoto-Nagai discs of mass respectively 10^" and 
10^^ Mq, of scale lengths respectively 1.5 and 4 kpc and 
common scale height of 0.5 kpc. The initial disc radius is 
25 kpc. 

Stellar velocities and velocity dispersions are found by 
solving numerically the Jeans equations directly on the 
grid using an iterative process. The initial conditions rep- 
resent an equilibrium solution of the stellar hydrodynam- 
ics equations, incorporating the softened gravitation as 
well as the natural anisotropic velocity dispersion tensor 
of disc. This obviously does not mean that the stellar disc 
is fully stable. Indeed, the disc is dynamically cold enough 
to be bar mode unstable. 

The gaseous component is represented by 50 000 par- 
ticles drawn from a third Miyamoto-Nagai disc of scale 
length of 6 kpc, scale height of 0.1 kpc and trunca- 
tion radius of 25 kpc. The total gas mass amounts to 
1.110^° Mq. An initial gradient of metallicity has been 
imposed following the relationship Z = 0.02 x exp(— r/r^) 
where r^, the scalelength of the metallicity distribution, is 
4.34 kpc. There is no vertical gradient as the gaseous com- 
ponent is initially distributed in a thin layer (scale height 
of 0.1 kpc). 

We also ran the same simulation with a factor of 5 
less particles to check whether the SFR and its variation 
in time depends on the number of particles. We find that 
our results reported in Sect. 13.21 are robust although the 
details of the dynamic evolution are obviously different. 
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For this first paper, we did not include any dark halo 
to be able in the near future to compare the effects of 
having added such an additional component. These initial 
conditions are identical to the Run A of Wozniak et al. 



Table 1. Calibration of the old population. 



After 1 Gyr, the simulation numbers 1 264 729 par- 
ticles, of which 26 573 gas particles. In comparison, the 
stellar disc and the dark matter halo of the final galaxy 
of Abadi et al. P003 ) is described by less than 166 000 
particles. The final galaxy of Samland & Gerhart (2003) 
contains 614500 stellar particles at the end of the simula- 
tion embedded in a growing halo. 

In Table 13 we give the masses of the various compo- 
nents (old and new populations, gas) for the initial time 
and after an evolution of 1 Gyr. 

2.3. Stellar mass-to-light ratio calibration 

The output of any N-body simulation is mainly a set of 
positions and velocities for all particles. In general, parti- 
cles have the same masses so that it is not necessary to 
make any assumption on their mass-to-light ratio. Thus, 
the problem of luminosity calibration has been widely 
overlooked. However, because of star formation, the par- 
ticles created during the run have various masses, ages 
and metallicities. The observable properties, which are 
mean quantities integrated along the line-of-sight (LOS) 
and over a spatial resolution element, must be weighted 
by the luminosity of the particles. Therefore, we need an 
accurate mass-to-light ratio to convert masses into lumi- 
nosities in various photometric bands. 

For each stellar particle, given its age and metallicity, 
the mass-to-light ratio T was obtained from a bi-linear 
interpolation into the tables of GISSEL2000 (Bruzual & 
Chariot ri993|l . for a Salpeter initial mass function (IMF), 
with mass cut-off at 0.1 and 100 Mq. Then, luminosity 
of the particle i (L^), in the waveband X, is computed 
according to its mass (M*), birth metallicity (Zl^^.^) and 
the time elapsed since its birth {t — t1^^^.^): 



^x{t ^born'^born) 



It is noteworthy that the photometric calibration is not 
fully consistent with star formation and feedback mod- 
elling since the selected value of i?"^^ corresponds to a 
lower mass limit for SNII equal to ~ 12 Mq whereas the 
standard value 8 Mq is used in GISSEL2000. 

For the old population (present at the beginning of 
the simulation) we have assumed ages of 6.7 and 10.4 Gyr 
and metallicities oi Z = 0.004 and Z = 0.02 (cf. Table Q). 
Obviously, this implicitly assumes that all the old stars 
simultaneously born 6.7 or 10.4 Gyr before the beginning 
of the simulation. Assuming an observer located at z = 
at the end of the simulation {t — I Gyr) , these calibrations 
imply a redshift of formation of z = 1 and z — S for resp. 
models C-D and models A-B. 

Global quantities (e.g. luminosities, mass-to-light ra- 
tios, etc.) and surface brightness distribution were com- 





Age 


Metallicity 


Model A 


10.4 


0.004 


Model B 


10.4 


0.020 


Model C 


6.7 


0.004 


Model D 


6.7 


0.020 



puted in various colours. We choose to work with four 
widely used broad band filters to cover the wavelength 
domain from the visible to the near infrared, namely B, 
R, H, and K in the Johnson's system of GISSEL2000. 

2.4. Dust extinction 

Dust extinction in B, R, H and K band is simulated in the 
simplest way, assuming a constant gas-to-dust ratio 

iV(Hi)/Av = 5.34 X lO^^cm"^ 

Ay is then converted to Ax where X stands for any wave- 
band, according to the standard interstellar extinction law 
from Rieke & Lebofsky H1985fl . H i column density should 
thus be integrated along the line-of-sight for each stellar 
particle, as Westera et al. (|2002|) did. Our SPH kernel is 
defined by (Monaghan & Lattanzio ll985|l : 



( 1-^ + ^ if0<f<l 



I 



if 1 < w < 2 



otherwise 



where v = r/h. The contribution of each particle to the 
density vanishes for r > 2h. Thus, for each stellar particle, 
we should have to draw up a list of gas particles of which 
the sphere of radius 2h is crossed by the line-of-sight. The 
gas column density is then the sum of individual gas par- 
ticle density weighted by the SPH kernel. 

This technique is obviously more difficult and CPU 
time consuming than with a grid. Therefore, we have com- 
puted extinction in a grid with the same spatial resolu- 
tion than our images and 11 slabs along the line-of-sight 
from z — —6.62 kpc to z = -1-6.62 kpc. The slab limits 
are spaced in log(|z|) in order to get a better resolution 
near the equatorial plane where most of the gas is concen- 
trated. Each slab absorbs the stellar luminosity behind it. 
For each slab, the gas density distribution is obtain by the 
convolution of particle positions by the kernel. 

3. Global evolution over 1 Gyr 

The simulation has been stopped at 1 Gyr for the following 
reasons: 

1. our simulation should be considered as a numerical 
experiment. Most of the physical assumptions and 
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recipes are simpler than in most chemodynamic nu- 
merical codes that are used to make simulations over 
several Gyr. This means that their validity might 
be questionable on longer timescales (e.g. the use of 
Toomre Q parameter). 

2. the closed box assumption is unlikely valid over 
a longer timescale. Gas infall, satellite harassment 
and/or merging, etc. are some examples of external 
perturbation that invalidate the closed box approxi- 
mation. 

3. since we do not add any dark halo (fixed or with live 
particles), we neglect the interactions between the disc 
and the halo. It is thus not reasonable to study the 
long term effect of the disc evolution without any dark 
matter. 

3.1. Dynamic evolution 

The initial disc quickly develops a typical strong bar and a 
spiral structure (Fig.P) both in the stellar and the gaseous 
components. The gravity torques due to the bar and spiral 
structure drive the gas inwards and reorganise the mass 
distribution even for the old stellar population; this gas 
inflow occurs in a rather short timescale since the star for- 
mation rate peaks att = 0.3 Gyr (cf. Fig.lSJ. As expected, 
star formation also occurs along gaseous spiral arms. 

After 1 Gyr, the total mass inside the central kpc has 
roughly doubled. Indeed, a nuclear gas disc is formed from 
the accumulation of gas in the centre, and new stars are 
actively formed there. The first cause of this mass inflow 
is the overall reorganisation of the mass distribution un- 
der the influence of the stellar bar. Due to the gravita- 
tional torques exerted on the gas by the stellar bar, the 
extra mass in the form of gas and new stars amounts 
to 3.5 10^ M0 at t = 1 Gyr, which is only 30% of the 
whole additional mass. Indeed, the old stars population 
contributes to the other 70%. 

The gas mass fraction inside the central kpc is very 
time dependent. For instance, between t — 0.905 Gyr and 
t = 1 Gyr the mass of gas is divided by a factor 4. This 
illustrates a property of the gas inflow towards the cen- 
tral kpc region: the fueling is not stationary but rather 
proceeds by burst. 

Kinematical properties of the central kpc are described 
by Wozniak et al. PMji . 

Table 2. Masses (in 10^° Mq) of the various population 
Sit t — and t ~ 1 Gyr. 



Population 


Whole simulation 
t = t=l 


Central kpc 
t = t = l 


Old stars 


11. 


1.04 


1.8 


Gas 


1.1 0.53 


0.08 


0.0011 


New stars 


0.57 





0.43 


Total mass 


12.1 


1.12 


2.23 




0.0 0.2 0.4 0.6 0.8 1.0 



Time (Gyr) 



Fig. 3. Global star formation rate (Mq yr""'^) versus the 
elapsed time from the beginning of the simulation. Note 
that star formation is inhibited during the first 60 Myr. 
Backward time average windows of 10 Myr (full line) and 
100 Myr (dashed line) has been used 



3.2. Star formation rate evolution 

The star formation rate (SFR) is displayed in Fig El The 
newly formed stellar mass is averaged backward over two 
time windows of widths respectively 10 Myr (SFRio) and 
100 Myr (SFRioo). Apart from the first SFRio peak at 
t = 0.1 Gyr which is due to the progressive switch-on of 
the star formation in the code, the other three SFRio max- 
ima occurring around 0.3, 0.5 and 0.6 Gyr are real. For 
i ^ 0.3 Gyr, the bar and a bi-symmetric spiral structure 
spontaneously form. The gravitational torques applied by 
the bar and the spiral arms on the gas create several re- 
gions of very high gas density in which star formation is 
ignited. SFRio reaches 30 Mq yr~^. As shown in previous 
works (e.g. Martin & Friedli [T997|l . the star formation is 
not homogeneously distributed over the whole disc but is 
mainly concentrated along the bar major axis and along 
the spiral arms (cf. Fig.j^J. Indeed, 3/4 of the star forma- 
tion occurs in the central kpc of the galaxy. SFR density 
ranges from 7.3 10"'* (in the disc) to 365 Mq yr~*kpc"^ 
(in the nucleus) . These values are in good agreement with 
data from Kcnnicutt H1998II . 

The next two peaks ai t k, 0.5 and 0.6 Gyr are the 
result of the gas inflow towards the central region of the 
disc. Such inflow is not continuous, but rather proceeds 
by burst. 

When the 100 Myr averaging window is used, the SFR 
maximum is moved forward by w 50 Myr. 
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Fig. 1. Snapshots of the evolution of the simulation from t = to i = 1 Gyr. From left to right are displayed the 
mass distribution and B dust-free, B with extinction, K dust-free, and K with extinction calibrated images. The field 
of view (20 kpc) is the same for each frame. The range in mass is [1., 1.3510^] Mq pc^^, displayed in log units. B 
surface brightness ranges from 12.75 to 28.8 mag arcsec"^ and K ranges from 10.9 to 26.2 mag arcsec~^. The particles 
have been rotated so that the bar position is roughly horizontal 



Fig. 2. Evolution of the simulation from t = to t = 1 Gyr. From left to right are displayed the gas mass distribution, 
the SFRio and SFRioo distributions and B— H colour map without and with dust extinction. The field of view (20 kpc) 
is the same for each frame. For the gaseous mass density, a lower cutoff has been apply for the sake of clarity. The 
range in mass density is thus [10~^, 2.5 10^] Mq pc~^, displayed in log units. SFRio and SFRioo ranges from 7.3 10""* 
(in the disc) to 365 Mq yr~^kpc~^ (in the nucleus), also displayed in log units. B— H ranges from —0.57 (black) to 
6.14 (white). The particles have been rotated so that the bar position is roughly horizontal 



3.3. Luminosity and colours evolution 

We plot in Fig. ^ the evolution of the integrated B— H 
colour index. Our values are in good agreement with 
Scodeggio et al. (2002^) . Indeed, for their sample of spi- 
rals, B— H values range between 1 and 4. For the sake of 
comparison, two SSP have been overplotted. One is the 
SSP of the old population and thus depends on the cali- 
bration model. The second is typical for a particle born at 
the maximum SFR (t = 0.3 Gyr). Its metallicity is Zq/2. 

In absence of dust, the galaxy becomes bluer by « 1 
magnitude between t = Q and t = 0.25 Gyr, then B— H 
slowly increases at a rate similar to the SSP of the popu- 
lation born at the maximum SSP. The colour of the disc is 
thus dominated by the population born during the most 
active phase of star formation around t = 0.3 Gyr. 

When dust extinction is taken into account, the am- 
plitude of the B— H evolution is strongly smoothed. B— H 
values remain close to the old population SSP. For models 
B and D, a B— H minimum can still be determined but it 



is shifted towards t w 0.4 — 0.5 Gyr. For model C, which 
has a young initial population with a low metallicity, B— H 
does not display any clear minimum. Indeed, the colour 
remains roughly constant over 0.6 Gyr. 

Fig- Ell shows the total luminosity evolution for all the 
models without (Ltot) and with dust extinction (L^^^). 
Independently of the old population calibration, Ltot evo- 
lution clearly shows two regimes. The broad peak centred 
around t = 0.3 Gyr obviously corresponds to the time 
when the star formation is the more active. There is a 
small but significant lag between the Ltot peak in the B 
and the infrared bands of « 10 Myr so that the B peak is 
correlated to the SFRio maximum while the K peak oc- 
curs between the SFRio and SFRioo maximums. The max- 
imum is followed by a slow decrease which corresponds to 
the dimming of the stellar population. However, the lu- 
minosity decrease is not exponential and is rather slower 
than for a SSP. This is due to a rate of star formation 
which remains high (above 2 Mq yr~^) during the next 
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Fig. 4. B— H colour evolution. Colours of the models with- 
out dust extinction are represented by full points, with ex- 
tinction by empty points. The full line stands for the B— H 
index of an SSP with the mean metallicity of gaseous star 
forming particles at t = 0.3 Gyr. The long-dashed curve 
is the B— H colour of the old population 

0.7 Gyr. This continuous star formation partly compen- 
sates for the ageing of the older stellar populations. 

Hereafter we call luminosity contrast the ratio of the 
maximum value of the luminosity to the luminosity at 
t = 0. For our calibration models, the latter corresponds 
to the luminosity of the old population at the beginning 
of the simulation, whereas the maximum value is the sum 
of luminosities of the old and of the new population. 

The Ltot contrast depends on the old population cal- 
ibration, especially in the B band. For instance, the Ltot 
contrasts are typically « 1.7 in K and between « 3.8 and 
« 5.6 in B according to the calibration model. However, 
these values are not typical for a SSP since the contrast in 
the infrared is always higher than in B. For instance, the 
contrasts amount to w 10 in K and « 5 in B for a SSP 
with Z = 0.004. This illustrates that the old population 
predominates in mass and that the B band luminosity of 
the old population is « 2 decade lower than the K band 
luminosity. 

As expected, when dust extinction is taken into ac- 
count, the total luminosities are fainter in all bands. 
However, since the amount of gas as well as its spatial 
distribution evolve because of the combined effects of the 
dynamic evolution and star formation, the ratio of un- 
obscured to obscured luminosity also evolves with time 
(Fig.Et). 

It is also noteworthy that the luminosity peak at 
t « 0.3 Gyr is smoothed out in R, H and K band. In B, 
LjQt still displays a broad bump with a small maximum 
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Fig. 5. a) Total luminosities for the four bands : B (blue), 
R (green), H (magenta) and K (red). Lines represent the 
total luminosity evolution of calibrated models without 
extinction (Ltot) while symbols stand for absorbed lumi- 
nosities (LjQt*^). Metallicity of the old population increases 
from left {Z = 0.004) to the right {Z = 0.02), its age 
decreases from top (10.4 Gyr) to bottom (6.7 Gyr). b) 
Evolution of the ratio L^^^ / Ltot in B (upper left panel), 
R (upper right), H (lower left) and K (lower right). The 
letters beside the lines refer to the calibration model 

but shifted after t > 0.4 Gyr. This delay is due to the 
fact that star formation with the highest SFR occurs in 
the regions with the highest gas density. In these regions, 
extinction is then the highest since we assume it is propor- 
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tional to the gas column density. Only between 1/4 and 
1/3 of Ltot is able to escape from the galaxy at i = 0.3 Gyr 
(cf. Fig.jSjs). Thus, the regions which contain the youngest 
population and thus contribute the most to the B band 
luminosity are also the most obscured (cf. Sect. I4.4.2|) . 
Since a fraction of the gas amount is consumed during star 
formation and since gas and stars have different kinemat- 
ics, the gas density decreases during the evolution. This 
enables the escape of the blue light after several tens of 
Myr. Consequently, the L^^^ maximum should not be in- 
terpreted as a SFR maximum since it occurs ~ 100 Myr 
later on. 

The high extinction of the youngest and brightest star 
forming regions also explains that the ratio Lj^^'^/Ltot is 
sensitive to the parameters of the calibration model (cf. 
Fig. ISJa) although the dust extinction model is exactly 
the same for all models. Indeed, since the B band lumi- 
nosity of young particles is absorbed by dust, the total 
B band luminosity is dominated by the contribution of 
the old population. The differences between Ltot and L^^^'' 
are more marked in B and R than in the infrared bands, 
since differences in luminosity are more pronounced in B 
and R band than in infrared bands for the different old 
population. But the mean extinction in K and H band is 
not neghgible since the ratio L^^^ /Ltot amounts to 0.7 at 
t = 0.3 Gyr. 

Stellar mass-to-light ratios are plotted in Fig. IHl 
Although the total stellar mass increases with respect to 
time because of star formation, the evolution of the total 
mass-to-light ratio is dominated by the total luminosity 
evolution. Values are in good agreement with mean values 
of the RC3 catalogue for spiral type galaxies (Roberts & 
Havnes [T994|l . Some observations (e.g. Vallejo et al. l2002|l 
present lower values for mass-to-light ratio in K band, in- 
dicating the presence of a very young population. 

3.4. Metallicity evolution 

Since an initial gradient of metallicity has been imposed 
on the gas particles at the beginning of the simulation 
(cf. Sect. I2.2jl . the metallicity of gas particles ranges from 
6 X 10"^ (at the disc edge) to 0.02 (at the centre). The 
mass-averaged metallicity (Z) is 7.5 x 10^^ (cf. Fig. d). 
Then, Z and the maximum metallicity increase until 
t = 0.3 Gyr. During the same time, Z of the new pop- 
ulation also increases and is always higher than for the 
gas since star formation is ignited in the central region of 
the disc. However, during the next 0.7 Gyr, the mass of 
the gas significantly decreases because of the strong star 
forming rate (SFRio~ 30 Mq yr~^). This implies that Z 
of the gas decreases while is remains constant for the new 
population. Some gas particles with a high metallicity are 
fully converted into new stars as their mass has reach the 
lower threshold. This explains that the maximum metal- 
licity of the gas decreases. 

At t — 1 Gyr, some low metallicity gaseous particles 
still exist, but because of the gas inflow towards the cen- 
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Fig. 6. Total mass-to-light ratio evolution. Lines, symbols 
and colours as for Fig. 
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Fig. 7. Metallicity evolution of gaseous particles (lines) 
and new population particles (symbols). Mass-averaged 
(continuous line and pluses), minimum (dot-dashed line 
and crosses) and maximum metallicities (dashed line and 
stars) are drawn 



tral region, their number decreases. When star formation 
becomes active for such particles (for instance when they 
encounter a spiral arm), they form new stars with a low 
metallicity, so that the minimum metallicity for the new 
population is a decreasing curve. 
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Fig. 8. B band surface brightness of one particle of the old 
population isolated in a cell of a CCD-like image at < = 0. 
The lines show the evolution of the surface brightness for 
the two metallicities used for calibration. The origin of the 
time is the beginning of the simulation (t = 0) so as the 
age of the old population is negative. The letters represent 
the model of calibration used for the old population (cf. 
Table nj. The 25'^ B band magnitude is plotted as an 
horizontal line for reference 

4. Surface brightness evolution 

4.1. Simulated images 

After individual particle calibration, simulated CCD im- 
ages are obtained summing particle luminosities into a 
512x512 pixels grid. The field of view is 60 kpc which 
gives a spatial resolution of ~ 117 pc, 1.3 times our small- 
est A^-body grid resolution. We thus produce one frame 
per waveband and per snapshot of the simulation. Our re- 
sults are obviously independent of the bar position-angle 
with respect to the North or any other axis. We thus de- 
cided to systematically rotate the positions of particles to 
align the bar with the x-axis. An example is displayed in 
Fig. n for model A. 

Due to a low number of particles in the external parts 
of the disc, we have also build two times lower resolution 
images. In Fig. |S1 we show the magnitude of the faintest 
pixel (in B mag arcsec"^) that is reachable in a simulated 
image, depending on the calibration model. This value 
corresponds to pixels which contains only one particle. 
Because this only happens in the outskirts of the disc, 
where star formation is inactive, these isolated particles 
belong to the old population. This figure also show that 
our mass resolution for the old population (2 x 10'°^ Mq 
per particle) disables the possibility to study calibration 
models with an initial age lower than 3.25 Gyr if we expect 



to determine an isophotal diameter at the 25* magnitude 
in the B band. 

To mimic real observations we should have to convolve 
our images with a point spread function. However this last 
stage depends on the telescope and site characteristics. It 
thus introduces a few free parameters which cannot be 
constrained without any detailed comparisons with real 
observations. 

In Fig. 121 B— H colour maps are displayed for model A, 
as well as the mass surface density of the gaseous compo- 
nent and the SFRio and SFRioo spatial distributions. In 
absence of dust extinction, B— H colour maps and SFRioo 
distribution are similar, which is no more the case when 
dust extinction is taken into account. The extinction is 
not uniformly distributed in the bar region. It is obvi- 
ously concentrated in the region of highest gas density 
and, therefore, in the region of high SFR. 



Each image has been analysed by fitting ellipses to the 
isophotes. This technique is very efficient in retrieving 
embedded structures (Wozniak et al. 119951 Friedli et al. 
I199fi|l . In general, it also gives reliable structural parame- 
ters (Wozniak & Pierce ni)91fl . The algorithm used in the 
ellipse-fitting program is fully described in Jedrzejewski 
(1987). Thus, description of the procedure will be brief. 

An isophotal contour of a given half-major axis length 
is written as: 

oc 

I{lp) = Co + c„ cos{n(p) + dn sm{nip) 

n=l 

If the trial ellipse exactly match the isophotal contour all 
the coefficients but cq vanish. Then, cq is the mean inten- 
sity along the contour. The departure from exact elliptic- 
ity can be quantified with the Fourier components: errors 
in the position of the centre of the ellipse give non zero 
values for ci and di while errors in ellipticity (defined by 
e = 1 — b/a where a and b are respectively the half-major 
and half-minor axis lengths) and/or position-angle (PA) 
are revealed by C2 and/or d2 components. 

Profiles of surface brightness, ellipticities and PA are 
obtained by increasing the half-major axis length a by a 
factor 1.01 between each fit. We have chosen this small 
value to obtain a good resolution in the inner region. 

The mean intensity profiles are then converted to sur- 
face brightness profiles. PA profiles follow the conventional 
notation. The North is at the top, the East at the left and 
the angles are counted in the anti-clockwise direction from 
0° to 180°. The bulk of the bar hes at PA« 0°. The range 
in radius for the figures begins at 1 px, i.e. 117 pc; the 
upper bound corresponds roughly to a surface brightness 
/iB ~ 10 mag arcsec^^. We show an example of these re- 
sults in Fig. ini for the model A at the same times than in 
Fig.d 



Time (Gyr) 



4.2. Ellipses fitting 
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Fig. 9. For each panel: from top to bottom, surface brightness (/i), ellipticity (e=l-b/a), position-angle (PA) and 
central shift for each fitted ellipse, for the two bands (B plotted in blue and H in red) and the mass distribution (in 
black). The mass density profiles have been scaled so that the 30*'* magnitude corresponds to 1 Mq pc^^. Results for 
dusty images are plotted in the form of dots 



4.3. Isophotal radii 

The surface brightness profiles obtained with ellipse fit- 
ting has been used to determined the isophotal radius in 
each band. The surface brightness levels at which isopho- 
tal radii are measured are 25 in B, 23.5 in R, 21.5 in K and 
20.5 in H. For the H band, we follow Gavazzi et al. ( 2000) 
while Gavazzi et al. (|1996|l used a limiting magnitude of 
21.5. 

Fig. EH displays the evolution of the isophotal radii for 
all the models. Apart from the H band, the isophotal radii 
increase with time during the first 0.6 Gyr and then re- 
main constant or slightly decrease until t — 1 Gyr. In the 
H band, R^q 5 clearly decreases after the SFR reaches its 
maximum value {t — 0.3 Gyr). The isophotal radius of a 
stellar disc, dynamically stable and free from star forma- 
tion, should decrease with time since at a given radius the 
surface brightness increases with time. However, the disc 
of our simulation shows a strong dynamic evolution. We 
have shown in Sect. 13. ll that the mass is fully redistributed 
under the influence of the bar. The mass inside the central 
kpc doubles after 1 Gyr and 70% of the additional mass 
come from the old population. Thus, the evolution of the 
isophotal radii could be the result of the combined effect of 
variations in luminosity due to star formation and stellar 
evolution, and the dynamic evolution. To disentangle the 
two effects, we plot in Fig.E]the evolution of the fraction 
of mass inside the various isophotal radii and, in Fig. 1121 
the evolution of the fraction of luminosity inside isophotal 
radii. 

The mass inside R^q 5 clearly decreases. But the total 
H band luminosity is rather constant or slowly increas- 
ing (Fig. Et) , while the fraction of luminosity inside R^q 5 
drops from 80 to 65%. Thus, the decrease of R^q 5 only 
reflects the dimming. It should be notice that R^q 5 is not 
a good estimator of the size of the galaxy since this radius 
only contains between 60 and 75% of the total mass and 
less than 80% of the total luminosity. The case of 5 is 
similar although the mass decreases slower. 



The evolution of the mass inside R25 and R23 5 is more 
tricky since it depends on the calibration of the old popula- 
tion. Rfg and R^j 5 are between 10 and 20% shorter when 
the metallicity of the old population is solar while they 
are only 10% shorter when the age of the old population 
is greater. These differences are roughly constant. Indeed, 
at these radii the old population dominates in mass and 
in luminosity. So that the differences in isophotal radii 
between models are directly linked to luminosity differ- 
ences between old population. Moreover these difference 
are roughly constant over 1 Gyr for a many gigayears old 
population. 

It is also obvious from Fig. ^| that dust extinction is 
negligible at these radii since the gas density is very low. 
The measurements of isophotal radii are thus robust. 

Finally, Fig. E| also shows that the isophotal radii are 
not a consistent estimator of the size of a galaxy. Indeed, 
depending on the waveband and the delay with respect to 
the maximum SFR, these radii delineate regions of various 
properties (luminosity, mass, etc.) which cannot be com- 
pared between galaxies and, for the same object, between 
different wavelength range. 

4.4. Bar characteristics 

4.4.1. bar length and isophotal radii 

As displayed in Fig.^and Fig. [3 a stellar bar forms, grows 
and becomes rounder during the simulation. Thanks to 
our calibration, we are able to measure the bar length at 
various wavelengths. Several criteria, defined in the liter- 
ature, can be used to determine the bar length. The re- 
sults can then be compared with dynamic radii. This work 
will be published in a separate paper (Michel-Dansac & 
Wozniak I2004|l . Hereafter, we will use bar lengths deter- 
mined with a criterion shown to be in good agreement 
with the corotation radius (namely C3 in Michel-Dansac 
& Wozniak I2004|l . The bar lengths being similar in H and 
K band, we only display the results obtained from the K 
band. 
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Fig. 10. Evolution of the isophotal radii: Rfg in B (blue), 
Rfs 5 in R (green), Rfo.5 in H (magenta) and R^i 5 in K 
(red). Symbols and lines as for Fig. |31 



Fig. 12. Evolution of the fraction of luminosity inte- 
grated inside the isophotal radii: Rfs in B (blue), R^3 5 
in R (green), R^p.s ^ (magenta) and R^i 5 in K (red). 
Symbols and lines as for Fig. [S] 
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Fig. 11. Evolution of the fraction of mass inside the 
isophotal radii: Rfg in B (blue), RI3.5 in R (green), R^q.s 
in H (magenta) and R^i.5 in K (red). Symbols and lines 
as for Fig. El 



In Fig. ^1 the evolution of the bar length to isopho- 
tal radius ratio is plotted. This ratio obviously increases 
during the bar formation phase then remains roughly con- 
stant and rather similar in both bands, for all the models 



apart from model B. Between t w 0.3 Gyr and i w 0.4 Gyr 
there is a peak. 

It should be notice that the difference between mod- 
els are mostly due to differences in the variation of the 
isophotal radius, since the length of the bar is almost in- 
dependent of the calibration model. 

The values are comparable to observational data from 
Martin p995|l . For a sample of isolated barred galaxies 
(Martinet & Friedli [T997|l . values of this ratio measured 
in B band range between 0.04 and 0.47. This ratio is in- 
sensitive to extinction effect since dust extinction does not 
corrupt measurements neither of the isophotal radii nor of 
the bar length. 

4.4.2. Luminosity of the bar region 

We have computed the total luminosity inside a circle of 
radius equal to the bar length. 

Comparison of our value with the literature is not 
straightforward. Indeed, our definition of the bar luminos- 
ity includes the bulge and the nucleus. Moreover, since the 
criterion used to determine the end of the bar is a good 
tracer of the corotation (cf. Michel-Dansac & Wozniak 
I2004f) . it gives upper limits of the bar lengths. In the case 
of real observations, the bulge and nucleus contributions 
are generally subtracted and bars end before the coro- 
tation. Thus, our values of the bar luminosity represent 
upper limits. 

Evolution of the bar luminosity to total luminosity ra- 
tio is plotted in Fig. E| for the four models, in B and 
K band. In B band, the ratio increases rapidly, reaches 
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Fig. 13. Evolution of the bar length to isophotal radius 
ratio: R25 in B (full circles), R21.5 in K (crosses) 



its maximum during the SFR peak, and then slowly de- 
creases. A few weak bumps correspond to secondary SFR 
peaks. According to the model, 80% or more of the total 
B band luminosity come from the bar region at the epoch 
of the peak of SFR. After 1 Gyr (SFR< 2M0 yr^^) more 
than two third of the total blue luminosity still come from 
the bar region. Extinction effects in the bar region are im- 
portant: at t w 0.3 Gyr less than 20% of the bar B band 
luminosity succeed to escape, and between 50 and 65% at 
i « 1 Gyr according to the calibration models. 

In the K band, the ratio increases rapidly, reaches 
a maximum (70-75%) during the SFR peak, and then 
rapidly (at t ~ 0.4 Gyr) reaches an approximatively con- 
stant value of about two third of the total infrared lu- 
minosity. Extinction effects are almost negligible, just as 
model effects. 

4.4.3. Bar ellipticity 

Maximum ellipticities reached in the bar region are dis- 
played in Fig. 1151 Bars are significantly thinner in B than 
in K band, except for t = 0.27 Gyr. In the B band, the 
bar is very thin until t w 0.55 Gyr. The maximum ellip- 
ticity then decreases with a slope which strongly depends 
on the calibration models. On the contrary, the K band 
maximum ellipticity decreases immediately after the SFR 
maximum. Hence, the difference between the two wave- 
bands increases until t « 0.6 Gyr, then it decreases with 
time but never vanishes. 

Two causes are simultaneously responsible of this ef- 
fect. Firstly, the spatial distribution of the new population 
is very elongated along the bar because new stars born in 



Fig. 14. Evolution of the bar luminosity to total luminos- 
ity ratio in B (full circles) and K (crosses). When extinc- 
tion is taken into account, the ratio is plotted as open 
circles in B and stars in K. 



the gas flow along the bar which is narrower than the 
stellar bar (cf. Fig. 12). Then, secular dynamic diffusion is 
responsible for progressively making it rounder. This ex- 
plains the high ellipticity reached in the bar when SFR is 
high and, only in part, the subsequent decrease. Secondly, 
the luminosity ratio between the new and the old popula- 
tions is wavelength dependent, begin higher in B than in 
K band. Moreover, it roughly follows the SFR variations. 
The colour dependence slowly becomes blurred as SFR de- 
creases. Hence, after the SFR peak, morphology becomes 
gradually dominated by the luminosity of the old popula- 
tion that has a rounder spatial distribution. This happens 
earlier in K band. In B band, the maximum ellipticity de- 
creases faster for more luminous old population (model 
C). This effect is weaker in K band and is linked to the 
age rather than metallicity of the old population. 

The maximum ellipticity of the mass distribution 
reaches its maximum at < = 0.27 Gyr, then slowly de- 
creases. The formation of the bar and its thickening oc- 
curs in the old stellar population which dominates the 
mass. However, the mass distribution is always rounder 
than the luminosity one. Since difference decreases with 
SFR, it could be anticipated that isophotes and mass iso- 
densities have the same shape if the SFR inside the bar 
vanishes. Thus, the surface brightness distribution, even 
in the near infrared, could be a good tracer of mass only 
if SFR remains low during a long time. Thus, mass mod- 
elling that used surface brightness distribution might give 
biased results if star formation is active, even at a moder- 
ate rate. 
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Fig. 15. Bar maximum ellipticity for all the models. 
Results for the mass distribution are displayed as black 
cross, red symbols for K band images and blue symbols 
for B band images. The letters beside the lines refer to the 
calibration model 



Dust extinction has obviously a more disturbing effect 
for B isophotes than K ones. Since extinction is propor- 
tional to H I column densities, it decreases with time fol- 
lowing gas consumption by star formation. Ellipticity pro- 
files are fully different and very noisy when extinction is 
high (cf. Fig.|5Jl, making maximum ellipticity determina- 
tion difficult. However, maximum ellipticity values do not 
significantly differ from the dust-free case^ . Differences on 
B band ellipticity profile are still present at t = 1 Gyr but 
appear mainly in the central kpc. 

5. Conclusions 

We have carried out a photometric calibration of a self- 
consistent N-body simulation of an isolated disc galaxy, in- 
cluding stars, gas and star formation. Dust extinction has 
been taken into account. We have simultaneously study 
the dynamic and photometric evolution in various wave- 
bands (B, R, H and K). We are able to reproduce obser- 
vational properties such as mass-to-light ratio and B— H 
colour indices. A morphological and photometric analysis 
have also been performed on simulated images in various 
wavebands, using the same tools as for real observations. 
Our main results can be summarised as follows: 



^ an error is slipped into Fig. 2 of Michel-Dansac & Wozniak 
ll^nn^ . On the left panel, we compare maximum ellipticity in B 
band with and without extinction, whereas on the right panel, 
we compare maximum eUipticity in B and H band. 



1. a dynamically self-consistent SFR peak produces a 
large increase in total luminosity in all wavebands 
whereas the stellar mass of the new population is much 
lower than the underlying old stellar mass. However, 
luminosity peaks can be almost completely obscured 
by dust extinction, since the young stellar popula- 
tion initially seats in gas-rich regions where extinction 
is maximum. This effect is particularly important on 
B— H colour maps since regions of highest SFR could 
become unobservable. It could lead to underestimate 
SFR from photometric or colour measurements. 

2. we show how isophotal radii evolve with respect to 
mass redistribution during the dynamic evolution but 
also with respect to stellar evolution. Hence, the dy- 
namic properties of the area enclosed by any isophotal 
radius depends on the waveband and on the star for- 
mation activity. Extinction effects do not affect the 
determination of isophotal radii, even in B band. It 
is also noteworthy that the luminosity integrated in- 
side the isophotal radius at 20.5 mag arcsec"'^ in the 
H band is not a good tracer of mass. 

3. the luminosity integrated into the bar length reaches 
« 80% of the total B band luminosity at the SFR 
maximum in absence of dust, and ^ 60% when SFR 
is low, whereas the bar length remains smaller than or 
equal to half the isophotal radius. 

4. when SFR is high (w 30 Mq yr'^) less than 20% of the 
B band luminosity of the bar region succeed to escape 
in the presence of dust because most star formation 
occurs in gas-rich regions where extinction is the most 
efficient. 

5. as long as star formation is active, mass isodensities 
of the bar region are rounder than isophotes. Thus, 
the surface brightness distribution, even in H band, 
is a good tracer of mass only when SFR is below 
1 Mq yr^^. Mass modelling that used surface bright- 
ness distribution might give biased results if star for- 
mation is active, even at a moderate rate. 

The use of other morphological classification tools (e.g. 
concentration and asymmetry, Abraham et al. will 
be will be reported in a forthcoming paper. 
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